MESOSCOPIC HYDRO-THERMODYNAMICS of PHONONS 



Aurea R. Vasconcellos (1) , A. R. B. de Castro (1 ' 3) , C. A. B. Silva (2) and Roberto Luzzi (1) . 
(1) Condensed Matter Physics Department-Institute of Physics "Gleb Wataghin", 
State University of Campinas- Unicamp, 
13083-859 Campinas, SP, Brazil. 
(2) Istituto Tecnoldgico de Aerondutica, Departamento de Fisica, 
12228-901 Sao Jose dos Campos, SP, Brazil. and 
(3) Laboratorio Nacional de Luz Sincrotron - LNLS, 13083-970 Campinas, SP, Brazil. 

Abstract 

A generalized Hydrodynamics, referred to as Mesoscopic Hydro-Thermodynamics, of phonons in 
semiconductors is presented. It involves the descriptions of the motion of the quasi-particle density 
and of the energy density. The hydrodynamic equations, which couple both types of movement via 
thermo-elastic processes, are derived starting with a generalized Peierls-Boltzmann kinetic equation 
obtained in the framework of a Non-Equilibrium Statistical Ensemble Formalism, providing such 
Mesoscopic Hydro-Thermodynamics. 

The case of a contraction in first order is worked out in detail. The associated Maxwell times are 
derived and discussed. The densities of quasi-particles and of energy are found to satisfy coupled 
Maxwell-Cattaneo-like (hyperbolic) equations. The analysis of thermo-elastic effects is done and 
applied to investigate thermal distortion in silicon mirrors under incidence of high intensity X-ray 
pulses in FEL facilities. The derivation of a generalized Guyer-Krumhansl equation governing the 
flux of heat and the associated thermal conductivity cofficient is also presented. 



I. INTRODUCTION 

It has been noticed Q that the ceaseless innovation in semiconductors design creates 
a demand for better understanding of the physical processes involved in the functioning 
of modern electronic devices, which operate under far-from-equilibrium conditions. The 
main interest is centered on the behavior of "hot" carriers (electrons and holes) but the 
case of "hot" phonons is also of interest, particularly in questions such as refrigeration of 
microprocessors and heat transport in small systems with constrained geometries 

These questions belong to the area of non-equilibrium phonon-dynamics jlj], more pre- 
cisely to the subject of phonon hydrodynamics associated to non-equilibrium (irreversible) 
thermodynamics [jfj], Q]. 

It may be noticed the fact that kinetic and hydrodynamics of fluids are intimately cou- 
pled. A first kinetic-hydro dynamic approach can be considered to be the so-called classical 
(or Onsagerian) hydrodynamics, which gives microscopic (mechano-statistical) foundations 
to, for example, the classical Fourier's and Fick's diffusion laws. But it works under quite 
restrictive conditions, (see for example js]) and, hence, more advanced approaches are re- 
quired to lift these restrictions, mainly, because, as noticed, the requirements for the analysis 
of situations involved in the nowadays advanced technologies. Several improved approaches 
were introduced as, for example, the Burnett approximation of hydrodynamics [9j, the one 
associated to Extended Irreversible Thermodynamics 6], and, recently, the so called Meso- 
scopic Hydro-Thermodynamics (MHT for short) dealing with heat transport [k| covering, 
in principle, all kinds of motion, that is, involving intermediate to short wavelengths and 
ultrafast motion. A complete MHT for classical fluids is reported in Ref. 

In this communication we describe the construction of a MHT of phonons, presenting an 
in depth study of phonon hydro-thermodynamics in the framework of a nonlinear quantum 



12|-[17||], referred 



kinetic theory based on a non-equilibrium statistical ensemble formalism 
to as NESEF for short. 

In section II derivation of a phonon higher-order hydrodynamics is presented. The 
hydrodynamic equations consist of the coupled sets of time-evolution equations for (a) the 
density of quasi-particles, together with those for its fluxes of all orders, and (b) the density 
of energy and its fluxes (heat transport) of all orders. The two families of equations are cou- 
pled together through thermo-elastic effects. The practical application of such encumbered 
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sets of equations requires to introduce a contracted description, that is, to retain only a 
finite number of fluxes and closing the set of equations by means of auxiliary complementary 
equations. 

In section III we fully derive a description said of order one, where we retain only the 
two densities and the first flux of each density. 

In section IV we present an analysis of some characteristics of these equations, namely, 



fa) the Maxwell times 



18j associated to the densities and fluxes, and (b) the construction of 



a hyperbolic Maxwell-Cattaneo-like equations for both densities and of a Guyer-Krumhansl- 
like equation for the heat flux. 

In section V we consider the simpler case when thermo-elastic effects are neglected and 
the equations for the two densities decouple. 

In section VI, on the contrary, thermo-elastic effects are analyzed and used to investigate 
transient thermal distortion of an optical substrate illuminated by a single intense ultra-short 
X-ray pulse as currently available in FEL facilities. 

Section VII contains our concluding remarks. Details of the calculation are given in 
several Appendices. 



II. PHONON MESOSCOPIC HYDRO-THERMODYNAMICS 

We consider a system of acoustic phonons in a semiconductor, in anharmonic interaction 
among them, in interaction with what we call a thermal bath of other degrees of freedom of 
the system, and in the presence of an external source capable of driving them out of thermal 
equilibrium. Moreover, the system is in contact with an external thermostat at temperature 
T . Since the divergence of a transverse field vanishes, and since in the absence of vortices 
the rotational also vanishes, in what regards the hydrodynamic motion we consider the 
longitudinal acoustic (LA) phonons only. 

The system Hamiltonian quantum mechanical operator is 

H = H os + H SB + H SP + Hqb (1) 

where Hqs is the Hamiltonian operator of the free LA system, namely, 

% = E q KK«q + l/2) (2) 
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with w q being the frequency dispersion relation, the wavevector q runs over the 1st Brillouin 
zone and a q , a q are annihilation and creation operators for LA phonons in mode q. 

The Hsb term accounts for the interaction of the LA phonons with what we have called 
the thermal bath, namely, the anharmonic interaction with the TA phonons, the effect 
of impurities elastic deformation and interaction with the electrons, or holes, they being, 
bonded or itinerant, the effect of imperfections (mainly dislocations), and we write for it 

Hsb = E q (A q i? q oJ + A q £+a q ) (3) 

where _R q and i? q indicate the corresponding operators associated with transition processes 
involving the interactions mentioned above; A q is the strength coupling in each case. For 
illustration let us consider the anharmonic interaction operator between the LA and TA 
phonons, namely 

Hsb = E qk {M qk a q 6j_ k &l k + M qk a q a^ k 6+ k + M q ' k a q a q _ k 6+ } 

+ Hermitian conjugate, (4) 

where the first contribution on the RHS describes the processes LA(q) — > TM(q') + TA(q"), 
the second the processes LA(q) -)> LA(q') + TA(q"), the third LA(q) + LA(q') -» TA(q") 
and the Hermitian conjugate describe the inverse processes. M qk are the appropriate ma- 
trix elements. The operators o q &q_ k &l k l e& d to the appearance, in the kinetic equation 
of evolution, of a linear term characterized by a relaxation time r* n , while the operators 
q "_ k 6^ k _and a q a q _ k & k produce non-linear terms sometimes referred to as Lifshits [20| and 
2 if ] contributions, respectively. It is worth mentioning that these non-linear terms 



a q a 



Frohlich 



produce, when the system is driven (by the external pumping source) sufficiently far from 



equilibrium, a so called "complex behavior" [22], 23J]]. In the case of phonons, it consists 



in the emergence of a type of Bose-Einstein condensation, propagation of long-lived solitons 



and a kind of "laser" action [[24J-[32|]. In what follows we shall take the case of not too in- 
tense excitation and disregard non-linear contributions, restricting the analysis to the linear 
(Onsagerian) regime, when the emergence of "complex behavior" is restrained according to 



33] 



Prigogine's theorem of minimum entropy production in local equilibrium 

Moreover, Hob is the Hamiltonian of the subsystems involved in the thermal bath, 
and Hsp is the interaction of the LA-phonons with the external pumping source, to be 
specified in each practical application. 



Next, to apply NESEF it is required first of a 
used to characterize the non-equilibrium ensemble [ 



1 to specify the basic dynamic variables 



14(- 12]]. A priori, when the system is 



driven away from equilibrium, it is necessary to include all observables of the system, which 



leads to the introduction of many-particle dynamic operators 
single-particle dynamic operator 

^qQ = «q~+Q/2^q-Q/2 



34|, [35]]. Here we use the 



(5) 



in the second quantization representation in reciprocal space, with q and Q running over 
the Brillouin zone. 

The two-phonon and higher-order dynamic operators can be ignored because of Bo- 
goliubov's principle of correlation weakening [ 36j, 37(]]. Moreover, since phonons are bosons, 
it would be necessary also to include the creation and annihilation operators a q and be- 
cause their eigenstates are the coherent states {38], and also the pair operators a^a^ (a^a^) 
because the number of quasi-particles is not fixed [39] . However, we disregard them because 
they are of no practical relevance for the problem at hands. In Appendix A we describe 
the corresponding non-equilibrium statistical operator. The energy of the thermal bath Hob 
is also a basic microdynamical variable and then we have the basic set 



(6) 



where, Q = corresponds to the occupation number operator z/ q = a^a^ describing a 
homogeneous phonon population and those with Q / 0, namely z/ q Q = «^ + Q/ 2 ^q-Q/2; 
account for changes in space of the non-equilibrium phonon distribution function. 

The average, over the non-equilibrium ensemble, of the microdynamical variables in 
the set of Eq. ([6]) provides the variables which characterize the non-equilibrium macroscopic 
state of the system. Let us call them 



{f q (t) =< Vq\t >; ^qQ(t) =< ^qQ|t >; E B =< H Q B >} 



(7) 



where Q ^ and < Vq\t >= Tr{d C{ Q £ (t)}, etc., that is, the average of the microdynamical 
operators of the set ([6]) over the non-equilibrium ensemble according to the formalism in 
Appendix A, where we have introduced the non-equilibrium thermodynamic state variables 
conjugated to those above, namely [cf. Eq. (lA3l) ] 



{F q (t), F qQ (t), A,}; 



(8) 



where l//3 — /csT . 

Going over to direct space we introduce the space and crystal momentum dependent 
distribution function 

^q( r > = ^Eq^qQ (t)exp[iQ ■ r], (9) 

V is the sample volume which shall be taken as one in what follows, where it may be noticed 
that f q (t), Q = 0, refers to the global properties and those with Q ^ reflect the changes in 
direct space. Hence, the hydrodynamic variables that describe the hydrodynamic movement 
consist of two families, namely 

{n(v,t), I„(r,t),...,/M(r,t),...} (10) 

with £ = 2,3,..., describing the movement of the quasi-particles, which we call the n-family, 
and where 

*M) = E q " q M) (ii) 

is the density of the quasi-particles, 

InM) =E q V q ^ q (r,t) (12) 

is the first (vectorial) flux of the density (current), 

^ ] (r,t)=E q V^ q ^ q (r,t) (13) 

is the rank-£, with £ > 1, (£ — 2, 3, 4 . . .) tensorial flux of the density, where 

w q = [V q w q VqWq...£ times... V q w q ] (14) 

meaning tensorial product of £-times the vector V q a; q , rendering a tensor of £-rank. 
On the other hand, we have 

{h(v,t), l h {v,t),...,lf{v,t),...} (15) 

again with i = 2, 3, . . ., describing the movement of energy of quasi-particles, which we call 
the h-family, and where 

Mr.^Eq^q^*). (16) 

Ifc(r, t) =Eq^ q V q Wq fq(r, t) , (17) 

4 £] (r,t)=E q ^qV^ q ^q(r,t), (18) 



which are, respectively, the density of energy, its first flux, and the higher order fluxes. 

By deriving on time t, the two sides of the general equations ffTTT) to (1131) and f|T6|) 
to f)18p . there follow, written in compact form, the hydrodynamic equations of motion 

|4 £] (r,t) = E q ^ ] (q)|- q (r,t) (19) 
where p is n or h (for the corresponding families), and 

^ 1 (q) = v^ q , (20) 

and, 

^ ] (q) = ^; q V^ q . (21) 

Evidently, all the hydrodynamic equations of motion are dependent on the equation 
of motion for the single quantity z/ q (r,t). Returning, for practical convenience, to reciprocal 
Q space, ^ q Q(t) satisfies the evolution equation 

J^ qQ (t) = Tr{(iti)- 1 [u (lQ ,H}g E (t)} (22) 

which is the average over the non-equilibrium ensemble (described by the statistical operator 
g e (t)), of the quantum- mechanical Heisenberg equation of motion for the microdynamical 
variable (operator) z/ q q given in Eq.Q. 

Direct calculation of the RHS in Eq. (1221) is extremely difficult and then it is neces sary 
to resort to the introduction of a more practical non-linear quantum kinetic theory [[3]- 17 ]] 
briefly described in Appendix B, which is applied using an approximation consisting in 
retaining only the collision integral which is second order in the interaction strength. The 
resulting evolution equation, when rewritten in direct r-space becomes a generalization of 
the Peierls-Boltzmann kinetic equation. As shown in Appendix B, once the limit of large 
wavelengths is taken, it acquires a form resembling the standard one, to be used consistently 
in what follows, namely 

d 

— z/ q (r, t) = V q w q ■ V r z/ q (r, t) + V q il q ■ V r z/ q (r, t) 

-r q ^ q (r,t)-^]+Z->,t) (23) 

v° = l/[exp(huJk B T )-l] (24) 



with 



where z/ q is the equilibrium LA phonon distribution at temperature To, Xq Xt ' describes the 
rate of change due to external sources/sinks acting on the system, T q is a reciprocal lifetime 
given in Eq. (125|) . II q is the self-energy correction of the LA phonon frequency giving uJ q of 
Eq.( l271) . which, if in Eq.([3]) we consider only the linear anharmonic interaction, that is, the 
first term on the right of Eq.(j4j) and its Hermitian conjugate, are given by 

T q = vr/^E k |M kq | 2 (1 + v£ A + 6 (fi k+q + a - c^ q ) , (25) 



TA 



l/[exp(m h /k B T )-l), (26) 



Wq = W q + Hq> ( 27 ) 

n q = (n/H 2 )P.V.J2 k |M kq | 2 (1 + u£ A + ^) / ( fik+q + Q k _ Wq ) ; (28) 

where fi k is the k-dependent frequency of TA phonons whose distribution is ^ k A , and P.V. 
stands for principal value. 

Inserting Eq. (l23l) in Eq. (119jl . the hydrodynamic evolution equations for the n and h 
families become 

§- t I?(r,t) = E q ^'(q) {V q a; q • V r z, q (r,t) + V q H q ■ V r z/ q (r,t) 

-r q [^ q (r,t)-^]+^ xt -(r,t)}. (29) 

Next, a closure for the set of Eq.( |29l) must be introduced. This means, that we 
must express the i/ q (r, t) appearing on the RHS in terms of the hydrodynamical variables 
belonging to the set. First thing to notice is that we are dealing with an enormous set of 
coupled integro-differential equations linking densities and fluxes of all orders. In reference 
j4o| (equation 7) an equivalent representation, in direct r-space, of these equations is given 
as 

I JM (r, t) + V r • 1^ (r, t) = Zjdh'C^ (r - r, t) if (r , t) + (r, t) , (30) 

where C * (r — r, t) if (r, t) stands for the contraction in all Z tensor indices. 

Equation ( 1301 is a continuity equation, having on the RHS a collision integral which 
accounts for sinks (system relaxation effects) and sources (pumping, driving the system out 



of equilibrium). It can be considered an extended Mori-Heisenberg-Langevin equation 41|. 
It represents a quite cumbersome set of coupled equations, of unmanageable pro por tion. To 
proceed further it is then necessary to introduce a contraction of description [42], that is 

8 



to say, a reduction in the number of basic quantities, retaining only a few fluxes. Hence, 
we must look, in each case, on how to find the best description using the smallest possible 
number of variables. In other words to introduce an appropriate contraction of description: 
This contraction implies in retaining the information considered as relevant for the problem 
in hands, and to disregard nonrelevant information [43]. 



Elsewhere 



42] we have considered the question of the contraction of description (reduction 



of dimensions of the nonequilibrium thermodynamic space of states), where a criterion for 
justifying the different levels of truncation is derived: It depends on the range of wavelengths 
and frequencies which are relevant for the characterization, in terms of normal modes, of 
the hydro-thermodynamic motion in the nonequilibrium open system. 

In other words, since MHT implies in describing the motion when governed by smaller 
and smaller wavelengths, or larger and larger wavenumbers, accompanied by higher and 
higher frequencies, in a qualitative manner we can say that, as a general "thumb rule," the 
criterion indicates that a more and more restricted contraction can be used when larger and 
larger are the prevalent wavelengths in the motion ( changes smoother and smoother in space 
and time). Therefore, in simpler words, when the motion becomes more and more smooth 
in space and time, the more reduced can be the dimension of the basic macrovariables space 
to be used for the description of the nonequilibrium thermodynamic state of the system. 

It can be conjectured a general criterion for performing contractions, namely, a con- 
traction of order r (meaning keeping the densities and their fluxes up to order r) can be 
introduced, once we can show that in the spectrum of wavelengths, which characterize the 
motion, predominate those larger than a "frontier" one, A? rr , +1 N = v 2 9 r 9 r+ i, where v is of 
the order of the thermal velocity and 8 r and 8 r+ i the corresponding Maxwell times 42]. 

In the next section we consider the situation when it is possible to use a contracted 
description including only the densities and their first fluxes. 



III. ANALYSIS OF THE PHONON MHT OF ORDER 1 

Let us consider a contracted description including the densities of quasi-particles n(r, t) 
and energy h(r,t) plus their first fluxes I n (r,t), I/i(r, t) only. For practical convenience it is 
better to work in reciprocal Q space. These quantities are then given by 

n(Q,*)=E^QQ(*), ( 31 ) 
9 



It is recalled that 



In(Q,0 =Eq V q^q 

h(Q,t) =E q ^ q V q u; q i/ qQ (t). 

qQ (*) = Tr { S q +Q/2 aq-Q/2 (*, 0) 



(32) 
(33) 
(34) 

(35) 



where (see Appendix A) 



g (t, 0) = exp{-cf>(t) - £ q [F q (i) £ q + E Q/ o^Q (*) ^qJ } (36) 
with, in this contracted description, 

^qQ (t) = (Q, t) + ifi h (Q, t) ^ + F n (Q, t) ■ V q W q + F h (Q, *) -^ q V q U; q . (37) 

Using the generalized Peierls-Boltzmann Eq. ( 1B6I) . the resulting equations of motion 
are (Q ^ 0) 

dn(Q,t)/dt = tQ ■ I„(Q,t)- (z/2) £ q (n q+Q / 2 - n q _ Q/2 ) z/ qQ (t) 

- (1/2) E q (r q+Q/2 + r q _ Q/2 ) Vqq (t) + xp xt - (Q, t) , 
dh(Q,t)/dt = zQ ■ i h (Q,t)- (z/2) E q (n q+Q/2 - n q _ Q/2 ) z/ qQ (t)^ q 

- (1/2) E q (r q+Q/ 2 + r q _ Q/2 ) z/ qQ (t)^ q + xf ext - (Q,t) , 

dI n (Q,t)/<9t = zQ-/P(Q,t)- (z/2) E q (n q+Q/2 - n q _ Q/2 ) u nQ (t)V 

- (1/2) E q (r q+Q /2 + r q _ Q/2 ) i/ qQ (t) V q u; q + xw ext - (Q, f) 



q w q 



di h (Q,t)/dt = zQ-4 2] (Q,t)- (i/2) E q (n q+Q/ 2 - n q _ Q/2 ) z/ qQ (t)^ q v 

- (1/2) E q (r q+Q/2 + r q _ Q/2 ) ^ qQ (t)fo; q V q u; q + zJJ ]ext - (Q, t) . 



q w q 



(38) 



(39) 



(40) 



(41) 



Next, we need to proceed to the closure of the equations, that is to express both 
z/ q q(t) and the 2 nd order fluxes in terms of the four variables n, h, l n and 1^. 

Regarding f q Q(£), we apply Heims-Jaynes [44j perturbative expansion for averages 
(around the homogeneous state) in a linear approximation, to obtain that, for Q ^ (see 
Appendix C) 

^qq(t) = h (q, t) n(Q, t) + b 2 (q, t) h(Q, t) + b 3 (q, t) ■ I n (Q, t) + b 4 (q, t) • I h (Q, t). (42) 
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Furthermore, developing n q _|-Q/ 2 around Q = and retaining only the lowest order 
non- vanishing term, that is, taking the long wavelength limit as indicated in Appendix B, 

(i/2) (n q+Q/2 - n q _ Q/2 ) « zQ-v q n q , (43) 
(i/2)(r q+Q/2 + r q _ Q/2 ) ^r q . (44) 

After introducing Eqs.fll2]), <$S§ and Eq.flEJ) in Eqs.flHEJ) to (JHJ) and, next, going 
over direct r-space we do have that 

dn(v,t)/dt + V-I n (r,t) =V-og (t) ■ I n (r,t)+V-afJ (t) ■ I h (r,t) 

+ 6n (t) n(r, t) + 6i 2 (t) fe(r, t) + Z^ ext - (r, f) , (45) 

dh(r,t)/at + V-I h (r,t) =V-a| (t) ■ I h (r,t)+V-a| (t) ■ I n (r,t) 

+ 6 22 (t) /i(r, *) + b 21 (t) n(r, t) + if ext ' (r, *) , (46) 

dlnfatydt + V-lMfat) =V-4 2 J (t)n(r,t)+V-a| (t) fc(r,t) 

+&g (t) ■ I n (r,t)+4 2 J (t) • I h (r,t)+Zf] ext ' (r, f) , (47) 

dl h (T,t)/dt + V-/f (r,t) =V-a| (t) /i(r,t)+V-4 2] (t)n(r,t) 

+&g (f) • I h (r,t)+&g (f) • I n (r,t)+Xl llext - (r, t) . (48) 

The coefficients afj (t), bfj (t), which depend on T q and on V q n q , are given in Appendix 
D. 

In order to close the system of equations, one still needs to express the 2 nd order 

[21 [21 

luxes In (r, t) and I h (r, t) in terms of the basic variables. We invoke again Heims-Jaynes 



44j perturbative procedure, in the linear approximation, to obtain that (see Appendix C) 

V.l|l(r,t) =B? ] (t) ■ Vn(r,t)+B? (t) ■ V/i(r,t), (49) 

V-i? (r,t) =cf (t) ■ Vn(r,t)+cf (t) ■ Vh(r,t), (50) 
Bf (t) = (q, t) [V qWq V q a; q ] , (51) 

Cf (t) = E q ^q&i (q, *) [V q w q V q u; q ] , (52) 

for j = 1 or 2, and coefficients (q, £) are given in Appendix D. Using the expressions 
above for V-In (r, t) and V-i| (r, i), the four equations, fl4"5]) to fl4"8"j) . become 

dn{r,t)/dt = [V-a'l (*) - v] •I n (r,t)+V-afJ (t) ■ I h (r,t) 
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+ b n (t) n(r, t) + b 12 (t) /i(r, t) + jM ext ' (r, i) 



(53) 



dh{r,t)/dt= V-4 2 J(t)-V -I^r^+V-a^ (t) • I n (r,t) 



dl n (r,t)/dt-- 



dl h (r,t)/dt-- 



+ b 22 (t) h(r, t) + b 21 (t) n(r, t) + xf ext ' (r, 










V-ag it) - Bf it) ■ V 


n(r,t)+ 


V-ag (t) - 


-4 2] 


it) 


V 


/i(r,t) 


+6| it) ■ I n (r,t)+4 2 J it) ■ I h (r,t)+ZW ext ' (r 


*)> 








v-41 it) - c? (*) ■ v 


/i(r,t)+ 


v-41 (*) - 




it) 


V 


n(r,t) 


+6S (*) " Ih(r,t)+6| (t) ■ I n (r,t)+4 1]ext - (r 


f). 









(54) 



(55) 



(56) 



The set of equations f )53|) to (1561) is a closed system of four linear first order differential 
equations for the four hydro-thermodynamic variables n(r, t), h{r,t), I n (r, i) and I^(r, i). 
Let us now analyze the contents of these equations. In Eg. (1531) (for the density n(r, t) of 
quasi-particles) and Eg. ( 1541 (for the density of energy h(r,t)), the first term on the RHS 
is the one of conservation, that is, the divergence of the corresponding flux; however, they 
are modified by the presence of contributions 41 (i) and 41 it) arising from the self-energy 
correction. The third term on the RHS corresponds to a relaxation-type contribution; the 
corresponding coefficients bn it) and b 22 (i) are minus the inverse of Maxwell characteristic 
times which we call 6 n it) and 9h it). They will be discussed in Section IV. The second and 
fourth terms are cross-contributions accounting for thermo-elastic effects coming both from 
phonon relaxation effects (coefficients b\ 2 it) and 621 it)) and phonon energy renormalization 
(coefficients 41 it) and 41 it))- The last terms Xn°' ext ' (r,t) and xj^ 6 **- (r,t) account for the 
rate of change generated by the externally applied driving agent. 

In Eg.( j55l (for the flux I„(r, t) of guasi-particles) and Eq. (l56]) (for the flux I/j(r, t) of 

[21 [21 

heat), the RHS terms with coefficient B- it) and C:- (£) had their origin in the divergence 
of the 2 nd order fluxes, which we expressed as superposition of the four basic hydrodynamic 
variables. The contributions with 41 it), 41 it), 41 (0 an d 41 it) are modifications arising 



J 2 ! U\ J 2 } 
u 42 

'■J 2 ] (+\ hi 2 ] 



out of self-energy corrections. The terms with 633 (i) and b A l it), similarly to the cases of 
n(r, t) and h{r,t), are relaxation-type contributions corresponding to minus the inverse of 
(tensorial) characteristic Maxwell times, which we call 8^ it) and 9^ it), see Section IV. 
As in the eguations for n(r,t) and h{r,t), we also find in the eguations for the fluxes cross- 
terms proportional to off-diagonal 's (due to self-energy corrections) and ft™ 's (due to 
relaxation effects), and external driving forces/sources. 
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Furthermore, we call the attention to the fact that in Eq. (l55|) for the flux l n (r,t) of 
quasi-particles, the first term on the RHS (roughly proportional to Vn) plays the role of a 
thermodynamic force analogous in classical fluid hydrodynamics to Fick's Law. In the same 
way, the first term on the RHS of Eq. (1561) for the heat flux I/ l (r,t) (the term here roughly 
proportional to Vh) is analogous to Fourier's Law. 

As a matter of fact, Eq. (l55l) and Eq. (l56l) . under steady-state conditions and after 
neglecting both external forces/sources and thermo-elastic cross-terms, reduce to 

In(r,t) = -D® (t) ■ Vn(r,t), where D® (t) = 

I h (r,t) = -Df (t) ■ Vh(r,t), where Df (t) = 

The .D's play the role of rank-2 tensor diffusion coefficients. 

In the next section we analyze several other aspects of this mesoscopic phonon hydro- 
thermodynamics of order 1. 



B?(t)-a£(t) -CW 



.PI 



3 [2] 



(57) 



Cf (t) - ag (t) -e®(t). (58) 



IV. CHARACTERISTIC MAXWELL TIMES AND MAXWELL-CATTANEO- 
LIKE HYPERBOLIC EQUATIONS 

We consider here some additional characteristics which can be derived from the 
results of the previous section. 



A. Characteristic Maxwell Times 



In Eqs. p5Tl - P8l the four coefficients b n , bf%, 622, &4J are minus the reciprocal of the 
so-called Maxwell times [47], namely 



- b n (t) = e- 1 (t) = E q r q &! (q, t) = V Wn (q, t) /r. 



- b§ (t) 



C 11 («> 



[2] 



E q r q [bs (q, t) V q c^ q ] = (q, t) /r q , 



- h 2 (t) = Ql 1 (t) = £ q r q & 2 (q, t) hw n = E q ^h (q, t) 



q- 



&S (*) 



0£ l (*) 



E„ F q t b 4 (q, t) V q W q ] fiw q = Eq^Ih (<1> *) / T < 



q- 



(59) 

(60) 
(61) 
(62) 



where T q is given in Eq. (l25l) . bj (q, t) for j = 1, 2, bj (q, t) for j — 3,4 and % (t) are given in 
Appendix D. l/r q = T q is the relaxation rate towards the equilibrium phonon population 



in mode q, which depends on q and u, 



q- 
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Eq.( )59|) to ( 162|) tell us that the Maxwell characteristic times are given by a 
Mathiessen-like rule involving all the relaxation times in each mode, r q , weighted by different 
kernels which are normalized, i.e., 

E q ^n (q, t) = E q w h (q, t) = 1, (63) 

E q «fi(q,*) = E q t«S(q.*) = i m , (64) 

where l' 2 ' is the unit diagonal tensor. 

As a consequence, if we assume all r q to be independent of q (not a possible physical 
situation, see below) then the Maxwell characteristic times are all equal. On the other hand, 
if we apply the mean- value-theorem of calculus, taking outside the E q a suitable mean-value 
of T in each case, say T^ h \ and r( Ih ), we find 6 n = r n = 1/r^, 6 h = r h = l/r( h ), 

efl = i^ Tln = ip]/r( ln ), 9 [ i = i [2 W lh = iM/r^). 

The quantity T q , given in Eq.(l25l). vanishes unless w q = f2k+ q + due to the 
5-function; then it can be rewritten identically as 



r q = E k l M k q | 2 5(n k+q + a 



U), 



{ (e^^ - 1) / (e^( q ) [1 - e -0M«O] + [l - e^^] } } (65) 

The matrix element |Mk q | 2 behaves as ~ |k + q|A;g. The frequency ui (q) of LA 
phonons vanishes at the center of the Brillouin zone and is maximum at the boundary of 
the zone. Then inspection of Eq.f l65|) allows us to estimate that T q is an increasing function 
of |q| over the Brillouin zone, therefore, r q = l/r q is a decreasing one with q increasing. 

The weighting functions in Eq.(l59l) to ( |62l) . and the time-dependent coefficients Aij 
and Ar. 1 below, are given in Appendix D, 

w n (q, t) = A7 1 (t) [A 22 (t) - Kw^A 12 (t)] z/ q (t) [1 + i/ q (t)] , (66) 



W In 



21 'q, t) = A3 1 (t) [V q u; q (4 2 J (t) - fc^gj (t))] • z/ q (t) [1 + z/ q (f)] V q w q , (67) 

w h (q, t) = A^ 1 (t) [A u (t) Hw^ - A 12 (t)] z/ q (t) [1 + z/ q (t)] hw^ (68) 

wfl (q, t) = A3 1 (t) [v q u; q (41 (t) hw^ - A| (t))] ■ ^ (t) [1 + i/ q (t)] /^ q V q o; q . (69) 
In a strictly Debye model it follows that 



# n (t) = e In (t) = £ q r> q (t) [1 + i/ q (t)] [(v 4 - 0Hs q v 3 ) 1 (p 2 p 4 - v\)] c-\ 
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(70) 



9 h (t) = 9 Ih (t) = £ q I> q (t) [1 + ^ (t)) [(-V 3 + (3hsqV 2 ) I (V 2 V 4 - Vf)] C-'PHsq, (71) 
where xd = /3HsqD, Qd is the Debye wavenumber and s is the sound velocity, 

V n = ! ° dx x n e x / (e x - l) 2 (72) 
Jo 

and 

C = (V/2tt) (q D /x D ) . (73) 

It must be stressed that the equalities 9 n = 6>i n and 9^ = ^ih are a consequence of using a 
strict Debye model (phonon group velocity independent of q). This can be a satisfactory 
approximation only under certain well-defined restrictions on the macroscopic state of the 
system. 

Summarizing, the characteristic Maxwell times 9 associated to the set of fluxes of all 
orders are composed by the weighted contributions of the relaxation times of the populations 
in each mode q, as described by Eq.( l59|) to ( 162|) . which are consistent with Mathiessen's rule. 
The 9 are all equal within a strict Debye model, but we stressed that this approximation is 
in general too restrictive. 

B. Maxwell-Cattaneo-like Hyperbolic Equations 

Let us consider the four equations, Eq.( |53|) to (1561) . where we recall that b u , b 22 , 
6 33 , 644 are minus the reciprocal of Maxwell times [cf. Eqs.f l59|) to f )62|) ]. and coefficients 

rol [0] 

Bf and Cf (j = 1,2) are given in Appendix D ( in a Debye model oj (q) = u (q) — sq, 
B [ 2 2] = Cf ] = and Bf ] = cf 1 = l®s 2 /3 ). After deriving in time Eqs.plD an d (El, 
next introducing Eqs.( l55l) and ( 1561) and assuming that the kinetic coefficients are weakly 
dependent on time, there follow the two coupled 2 nd order differential hyperbolic Maxwell- 
Cattaneo-like equations: 

d 2 n(r, t)/dt 2 + (d£ + 0" 1 ) dn(r, t)/dt + (- V-Sf 1 (t) • V + 9^9 n 1 + 634&21) n(r, t) 

+ (& 3 4 + 621) dh(r, t)/dt + (- V-Bf (t) • V + O12 + O34) h(r, t) 

=S n (r,t), (74) 
d 2 h(r, t)/dt 2 + (6£ + 9 h x ) dh(r, t)/dt + {-V-Cf it) ■ V + 9^9 h 1 + 643&12) h(r, t) 
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+ (&4 3 + 621) dn(r, t)/dt + (-V-Cf 1 (t) • V + O21 + O43) n(r, t) 

=S h (r,t). (75) 

In these equations, the contributions associated with self-energy corrections were neglected 
and and account for sources and/or external forces. 

V. DECOUPLED MOTIONS OF QUASI-PARTICLES AND HEAT 



Consider again the set of Eq.( l53l) to (1561) . in real r space, for the four thermodynamic 
variables n, h, I n and 1^. From now on we neglect the coefficients a^, i.e., we disregard contri- 
butions arising out of self-energy corrections, see Appendix D. We also neglect coefficients 
bij with i ^ j. Such describe thermo-elastic effects, which will be discussed in Section 
VI. 

The set of equations Eq. (153]) to (155]) becomes simplified and quantities related to 
the family n become decoupled from the quantities related to the family h and vice-versa, 
giving the two pairs of equations: 



dn(r,t)/dt+V-I n {r,t) = -9^ (t) n(r, t) + Xp xt - (r,t) 



dl n (r, t)/dt+B? ] (t) ■ Vn(r,t) = - [0£ (t)] 121 • I n (r,t)+jW ext - (r, t) 



-1 Ml PI 



(76) 
(77) 



dh{r, t)/dt + V.I A (r,t) = - Q- h x it) h(r, t) + Xf e (r, t) , (78) 

dl h (r, t)/dt+C® (t) ■ Vh(r,t) = - K (*)] m ■ Ih(r,t)+Xl 1]ext - (r, t) . (79) 

The coefficients present in these equations depend on time but not on position r. Then, tak- 
ing the time-derivative of Eq. (l76|) and the divergence of Eq. (l77j) one can partially eliminate 
the flux I n ; the intermediate equation for n(r,t) is 



d 2 n(r, t)/dt 2 - V • \b [2] (t) ■ Vn(r,*)l = -d [O' 1 (t) n(r, t)] /dt + «9X™ ext - (r, t) /dt 



+ V 



C (t)l [2l -I n (r,t)-V-XW ext - (r,t). 



(80) 



If, in addition, the tensor 



l[2] 



is isotropic 



[2] 



0in (*) = C J (*)1 [2] )> the flux I n 



can be completely eliminated with help of Eq. (l751) . giving 

d 2 n(r, t)/dt 2 - V • \B [2] (t) ■ Vn(r,t)] = -d (t) n(r, t)] /dt + «9X,M ext - (r, t) /dt 
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+ t 11 (*) {-M*, t)/dt - 6- 1 (t) n(r, t) + X™ ext ' (r, t)} - V-jW ext ' (r, t) . (81) 
If (t) is isotropic (B^ (t) = Bil^) , the equation becomes even simpler 

d 2 n{r, t)/dt 2 - B 1 V 2 n{r,t) = - d [9~ l (t) n(r, t)] jdt - ^ n 1] (t) dn(r, t)/<9t 

+ C 11 (0 ^ ]ext - (r, f ) + dl^- (r, f ) /0f - V-XW ext - (r, f) . (82) 

The physical interpretation of the terms in this hyperbolic 2 nd order partial differen- 
tial equation for n(r, t) is as follows: The LHS has the form of a standard wave-equation in 3 
dimensions; the square propagation speed is B\. On the RHS we find two damping terms (the 
ones with time derivatives of n(r, t)) which depend on the two characteristic Maxwell times 9 n 
and #i n , (see Section IV) and a sum of three other terms which depend on derivatives of the 
given sources/external forces, namely^" 1 ' (t)X^ ext ' (r, t)+<9li°' ext ' (r,t) jdt— V-X^ 6 **' (r, t). 
The sum of these three terms is the effective source in Eq. (1821) . 

With exactly the same procedures, we get for h(r,t) an equation of the same form, 
but where all n are replaced by h: 

d 2 h(r, t)/3t 2 - C 2 V 2 h(r,t) =-d [9 h x (t) h{r, t)] jdt - 9 [ Ih 1] (t) dh{r, t)/dt 

+ 4 1] (*) lh ]eXt ' (r, *) + <9X? ext - (r, t) jdt - V-Xi 11 ^' (r, t) . (83) 

It can be noticed that these hyperbolic equations resemble the so-called telegraphist 
equation in electrodynamics. Regarding the solution of Eq. fl82l) and flHBl we can state that 
whenever the effective source has no spectral components at the frequencies of the hydro- 
dynamic modes, there is a unique "particular" solution obtainable with a Green's function. 



Morse and Feshbach 45] give such Green function for the simpler case of a strongly localized 
effective source. However, if the effective source has spectral components at the frequen- 
cies of the hydrodynamic modes, "particular" solutions may still be found but they will in 
general diverge as \t\ — > oo. 

Although the equation for n(r, t), Eq. fl82|) . and the equation for h(r, t), Eq. fl83l) . have 
the same form, the coefficients may be widely different in practice. For instance, in Eq.( l82l) 
for the phonon density n(r, t), the damping (terms with dn(r,t)/dt) is related to sound 
attenuation in the material medium, which is very small in most rigid material media. In 
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this case, the damping terms are a small correction on the undamped wave equation and 
can frequently be ignored. 

Regarding Eq. ( 1551) for the density of energy, suppose conditions of quasi-equilibrium 
have been attained. Then we can define a local quasi-temperature T*(r, t) such that h(r, t) = 
CyT*(r,t) where Cy is a constant specific heat per unit volume. Then, for steady-state 
processes in ordinary solid material media, it is empirically established that d 2 h(r,t)/dt 2 is 
negligible in comparison with the damping terms dh(r,t)/dt. In this case, Eq. (l83l) reduces 



to Fourier's equation for T*(r, t) 46]. 



Let us consider heat motion in more detail. It is governed by Eq.( l78l) and ( 1791) . 
Deriving Eq. (1791) in time and taking the spatial gradient of (178 p . we eliminate h(r,t) and 
get the following hyperbolic equation for the heat flux I^(r,t): 

d 2 l h (r, t)/dt 2 + + dl h (r, t)/dt + e^I h (r, t)+C 2 V (V-I,(r, *)) 

= a4 1]ext - (r, t) /dt + C 2 VjJ° ]ext - (r, t) . (84) 

We recall that V (V-I h (r, t)) = V 2 I /i (r, t)+V x V x I h (r,t). In that way, Eq.([84D is an 
extended version of the Guyer-Krumhansl equation, which in the steady-state and assuming 
V x I/i(r, t) = reads as 

I,(r)+^V 2 I,(r) =^VX? ext - (r) , (85) 

where, 



c 2 e h e m , (86) 



with having dimension of length. Notice that in the Debye model, C 2 is one third of the 
square of the sound velocity. 



VI. THERMO-ELASTIC EFFECT 



When a pulse of energy, well localized in space and time, is absorbed by a solid, 
the solid will heat and thermally expand, creating local time-dependent stresses and strains. 
As time elapses, the absorbed energy will spread away from the "hot spot" and the solid 
will eventually reach a state of thermal equilibrium with a uniform temperature, as well as 
mechanical equilibrium with vanishing stresses and strains. 
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Of immediate practical relevance in this analysis are the scalar field T*(r,t) describ- 
ing a local quasi-temperature, and the vector field u(r, t) describing time -and position- 
dependent displacements of the atomic nuclei. 

Consider a homogeneous equilibrium solid with uniform density £ - If now the solid is 
locally excited in some way and the nuclei in some neighborhood undergo time dependent 
displacements, the density, £ (r,t), will also become time -and position- dependent and for 
small displacements can be written as |47j 

£(r,*) = £ [l-V-u(r,t)]. (87) 

It is clear that all motions of the nuclei can be expressed in terms of the normal modes 
of vibration in the solid, or, equivalently, in terms of the phonon density, be it crystalline or 
vitreous. Hence we argue that, upon a time -and space- localized excitation in a solid, the 
change in material density is, within a multiplicative constant, the same as the change in 
n(r,t), the phonon-density discussed in the previous sections. Therefore, we write 

Vn(r,t) = -£ V[V-u(r,t)], (88) 

dn(r, t)/dt = -£ V ■ <9u(r, t)/dt, (89) 

I n (r,t) =£ du(v,t)/dt. (90) 

If we neglect the off-diagonal coefficients ay and fey in Eq. fH5|) and Eq. (|47|) and use 
the definitions of Maxwell characteristic times (see Eq. (1591) and Eq. (l60l ). the former become 

dn(v,t)/dt = -V-InM-C 1 (i)n(r,t) +Z™ ext - (r,t) , (91) 

dl n (v,t)/dt = -V-I®(r,t) - e£ {t) I n (r,t)+jW ext - (r,t) . (92) 

The five equations (188]) through (192]) can be combined to produce the hyperbolic 
equation 

£ d 2 u(r,t)/dt 2 + fa/fa) d*(T,t)/dt = -V-/M(r,t)+XW ext - (r,t) . (93) 

To close this equation, we need to express V-iV (r, t) in terms of a mixed representation, 
what is described in Appendix C. We end up with 

d 2 u(v,t) „ i . . du(r, t) _ -i . r9l , . „ r „ 

d \ 2 + 6 ln (*) + ^ A[ nn (*) ' V [V • U (r, t)\ = 
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= -#W4J (*) • VT* (r, t) + to l $ ]ext - (r, t) . (94) 
Moreover, as shown in Appendix C, the evolution equation for the nonequilibrium tem- 
perature (quasitemperature T* (r, t)) is 

e * d 2 T*(r,t) | dT*(r,t) | 0* fc(r,t) 



«9t 2 <9t 9 h 9 lh C v (r,t) 



B n — — B h 



8* 1 Xr Xt '(r,t) 



Where, 



* B T. 2 C,(r, ( ) Vr(r - () + ^T- (95) 
+ (96) 

^ = Eq (^q) 2 Vq^q ■ Vq^q 5Vj5F n , (98) 

with I7 q (t) given in Eq. flClip . 

In the evolution equation for u(r, t) Eq.f l9~4j) we find in the RHS a term which is 
)roportional to VT* (r, t) in agreement with the phenomenological equation of elasticity 



471 ]. In the application to be considered here, there is no external source of "particle flux" 
and we set Z^ ext ' (r, t) = 0. In addition, the term du(r,t)/dt describes sound attenuation, 
which will be negligible in the illustration, hence we drop it. 

In the evolution equation for T* (r, t), Eq.( l95i) . there is an external source feeding 
energy in the solid medium, namely, a radiant energy pulse, so we keep the term xj^ ext ' (r, t) 
to be further specified below. Except at very short delays after excitation, the flow of heat 
is diffusive, which means that d 2 h(r,t)/dt 2 < {d h l + 6^) dh(r, t)/dt, therefore we drop 
the 2 nd order time-derivative. The effect of the 2 nd order time-derivative at short delay 



times was investigated in 50|. The term linear in h(r,t) describes time-relaxation effects 
which we lump together with the effect of the I s * order time-derivative. We also neglect 
the term with n(r,t), because n is almost uniform and heat transport is weakly affected by 
small fluctuations in material density; in other words, the measured macroscopic diffusion 
coefficient Dh already contains in itself the effect of density fluctuations which will be always 
present. 

Given these considerations, we take as practical equations the following: 

dT* (r, t) fdt - D h V 2 T* (r, t) = Cy l w{r, t), (99) 
20 



d 2 u(r, t)/8t 2 - c 2 s V (V ■ u(r, £)) = - (a v K/£ ) VT* (r, t) , (100) 

where D h is the diffusion coefficient, Cy is the specific heat per unit volume, w(r, t) is the 
power density being transferred from the radiant pulse into the medium, c s is the speed of 
sound, ay = (l/V)(dV/dT) is the volumetric coefficient of thermal expansion, K is the bulk 
modulus and £o is the material density. 

Now we can address, as illustrative application, the transient thermal distortion of 
an optical substrate (mirror, Bragg crystal, diffraction grating, etc..) illuminated by an 
intense ultra-short X-ray pulse such as currently available at several Free-Electron-Laser 
(FEL) facilities in the world jsi ] . This situation is entirely different from the case of optical 
components subjected to steady-state heat loads. 

The pulses produced at FEL facilities are strongly localized in 3-dim space and in 
time. When such a pulse of X-rays reaches the interface vacuum/solid at the surface of an 
optical element, some fraction of the radiant energy is reflected and some fraction penetrates 
the solid to a depth So and is absorbed, generating heat and thermal distortion which impacts 
the optical performance of the device. We want to estimate the seriousness of the surface 
distortion in a time scale of pico- to nano-seconds after incidence of a single X-ray pulse 
lasting only a few femto-seconds. The analysis below will show that, although the incident 
femto-second pulse is gone long before the optical surface has time to distort, the next pulses 
in a pulse-train can be badly affected. 

Consider the interface between vacuum at z < and an infinite slab of silicon at 
z > 0, extending indefinitely in the x and y directions. Crystalline silicon is a very popular 
material for optical substrates, because it is available in large sizes, it has high thermal con- 
ductance and low coefficient of thermal expansion (hence thermal distortion is minimized), 
accepts state-of-the-art polishing (RMS surface roughness of only a few Angstroms) and is 
relatively cheap. 

A short Gaussian-shaped pulse of X-rays (duration t FEL , radius r ) centered at 
wavelength Xfel and propagating along the z axis comes from z = — oo and hits the 
interface at t — 0. This scenario is cylindrically symmetric about the z axis. Tables I and 
II describe the FEL pulse and the solid medium. 
Table 1 - FEL characteristics 
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Quantity 


Symbol 


Value 


Unit 


Photon energy 


hi> = fvjj 


5000 


eV 


Photon wavelength 


A 


22.8 


A 


Pulse duration 


/\ T, T? 77 T 


10 


f sec 


Pulse total enerEfv 


Wtptpt 


0.8 


TnJ 


Peak FEL power 


IIP TP T 


80 


GW 


4^ nhotons ner FEL nulse 




1.0 x 10 12 




Optical reflectivity 


R 


0.8 




Absorbed energy 


(1 - R)W FEL 


0.16 


mJ 


Pulse Gaussian radius 




2.0 x 10 6 


nm 






2.0 


mm 



Table 2 - Thermo-elastic constants for silicon 



Quantity 


Symbol 


Value 


Unit 


Density 




2.329 


g/cm 3 


Specific heat (per unit mass) 




0.702 


J/(g °c) 


Specific heat (per unit volume) 




1.635 


J /(cm 3 °C) 


Thermal Conductivity 




1.68 


W/ {cm °C) 


Heat diffusion coefficient 


D h 


0.102 


nm 2 / f sec 






1.02 


cm 2 /sec 


Coeff. of ther. exp. (linear) 


Oil 


3 x 10" 6 


1/°C 


Coeff. of ther. exp. (volumetric) 


a v = 3ai 


9 x 10~ 6 


1/°C 


Elastic (bulk) modulus 


K 


1.06 x 10 12 


g/ (cm sec 2 ) 


Shear modulus 


c 


2.19 x 10 11 


g/ (cm sec 2 ) 


Poisson's ratio 


V 


0.45 




Speed of sound 


C s 


7.21 x 10~ 3 


nm/ f sec 






7.21 x 10 5 


cm/ sec 


Opt. abs. const, (hu = 5000 eV) 




0.555 x 10" 4 


l/nm 



The power density w(r,t) to be used in Eq. (l99|) is 



w(r,t) = 0, (101) 
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if z < 0, for any time — oo < t < oo, 



w(r, t) — P (x, y, z, t) exp (— fjiz) , (102) 
if z > 0, for any time — oo < t < oo, 

P (x, y, z, t) = Uon 1 / 2 (z /t FEL ) exp [- (x 2 + y 2 ) /r 2 ] 5 (z - ct) . (103) 

For simplicity we take the limit of "short" z and tpEL, with (^oAfel) = c, the 
speed of light in vacuum. Integrating P over all time and all space we get the total energy 
Wfel = Uq tt^VqZo carried by the pulse, and find that Uq has the meaning of average 
energy density in the radiant pulse, /i is the light absorption coefficient in the solid medium, 
which, for a given medium, depends on the wavelength Xfel of the radiation. Notice that 
the calculation is linear and all results scale linearly with the amount of absorbed energy, 
which in this application is 0.16 mJ. 

The procedure here will be the following. First, one solves Eq. (l95|) for the tempera- 
ture field T*(r,t). Next, one uses VT* (r, t) as source in Eq. flMj) for u(r,t). 

Green's functions, analytical expressions for the solutions T*(r, t) and u(r, t), and 
a discussion of numerical methods are given in Ref. j48|. Here, let us recall that u(r, t) = 
u Par *(r, t)+u Free (r, t) where u Pari is a particular solution which depends on the given source, 
while u Free is any arbitrarily chosen solution of the associated homogeneous equation (no 
sources), chosen to satisfy boundary/initial/asymptotic conditions, as the case may be. The 
condition to be met here is that all normal stresses at the "free" surface z = of the solid 
medium be vanishing: a zi (z = 0) = a iz [z = 0) = 0. These conditions lead to coupled 
Fredholm integral equations of I s * kind, which we have solved only approximately using 



truncation, see Ref. [481 ]. 



The result of these thermo-elastic calculations is that, on absorption of 0.16 mJ of 
energy from a Gaussian X-ray photon pulse at Huj = 5000 eV, with r = 2.0 mm, lasting 10 
fsec, there is an outwards surface bulge several nm high that comes after the FEL pulse, 
with a delay of several hundred nanoseconds. 

Detailed results are shown in Figures 1 to 5. 

Figure 1 shows the temperature T* versus depth z inside the Silicon slab, at selected 
times after incidence of the X-ray pulse, assuming T* = °C as initial temperature. The 
thick vertical line is the "causal cut-off" at t — 100 /sec; namely, for t = 100 fsec, the light 
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FIG. 1: Figure-1 Temperature T* versus depth z inside the solid medium, at p = 0, at selected 
time-delays t after incidence of the pulse, assuming T*(r,t < 0) = for all r. 



FIG. 2: Figure-2 Surface displacement u^ art {r = 0,t) versus time-delay t after incidence of the 
pulse. 

penetrates only as far as ct = 2.998 x 10 4 nm, hence the temperature at \z\ > ct is still 
identically zero. The cutoff for the other is off-range in this figure. However, as t increases, 
the X-ray pulse penetrates deeper and deeper till it is depleted by absorption. Then, further 
heating of Silicon layers far away from the surface proceeds by diffusion only. 

Figure 2 shows the surface displacement of the Silicon slab, at the center p = of 
the X-ray light spot, as a function of elapsed time. This is a contribution of the particular 
solution. Negative displacement means a surface bulge. The bulge is maximum at about 
400 nsec, and its amplitude is very large about 6 nm then goes back to the original position 
on a time scale of msec. 

Figure 3 shows the radial dependence of the surface bulge at selected moments. From 
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FIG. 3: Figure-3 Radial profile of the particular solution u^ art (p,z = 0,t) at selected time-delays 
t after incidence of the pulse (logarithmic vertical scale) 



FIG. 4: Figure-4 Radial profile of the particular solution u^ art (p, z = 0, t = 400 ns ) (linear vertical 
scale), for comparison with Figure 5. 

the profile at 400 nsec we find a surface figure error (maximum slope of the surface) of about 
2 [irad, which is not small, and takes microseconds to decay. Notice that the first X-ray 
pulse hitting "cold" silicon surface suffers no adverse effects, because at the time the surface 
bulges out, the pulse is already long gone. If, however, the experiment envisages a train of 
FEL pulses, and the spacing in the train is a few hundred nsec, the pulses following the 
first will be defocused by the heat-induced surface bulge. Furthermore, the bulge can be 
resonantly enhanced in a disastrous way. 

Figure 4 shows again, for t = 400 ns, the radial profile of the surface displacement 
u Part p rec ij c ted by the particular solution but in a linear scale for easy comparison with 
Figure 5. 

Figure 5 shows, for t = 400 ns, the radial profile of the surface displacement uf ree 
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FIG. 5: Figure-5 Radial profile of the free solution u Free (p, z = 0, t = 400 ns) (linear vertical scale), 
for comparison with Figure 4. The range < p < 7 mm was cut into several overlapping segments 
indicated by the differente symbols. The lack of perfect overlap at the edges of neighbouring 
segments gives an estimate of the errors incurred in the approximate numerical solution of the 
coupled integral equations. 

predicted by the "free" solution. This was obtained by numerical methods, after truncation 
of the coupled Fredholm integral equations of I s * kind which follow from the requirement 
of vanishing normal stresses at the surface z = 0. This has opposite sign and is smaller 
than the particular solution. The range < p < 7 mm was cut into several overlapping 
segments in order to speed up convergence; different symbols indicate distinct segments. 
The lack of perfect overlap at the edges of neighbouring segments gives an estimate of the 
errors incurred in the approximate numerical solution of the coupled integral equations. 
The net displacement is wf art + u F / ee = — 5 nm (the minus sign means the displacement is 
outwards). The slope du/dp is known as "surface figure error" and is an important figure 
of merit for optical components. Here it is ~ 2 rad, which significantly exceeds the 
current state-of-the-art in optical polishing. 

We can compare the present results to the previously studied case of vacuum ultra- 



violet light [48|], |49|], [50J where the photon energy was ftuj = 100 eV, absorbed energy 
per unit area Wq = 180 piJ/cm 2 , absorption coefficient 0.128 (nm) -1 , peak displacement 
UMax = 0.02 nm and "specific displacement" UMax/W = 1.3 x 10~ 4 nm/ '(/iJ / 'cm 2 ). 

For X-rays, present calculation, we find UMax/Wo = 8.4 x 10~ 3 nm/(p,J/cm 2 ), which 
is 65 times larger, even though the absorbed energy is comparable. However, if /i is the 
wavelength-dependent absorption coefficient, the penetration depth is of order l//i. This 
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is the thickness of the layer that is strongly driven out of equilibrium by absorption of the 
radiant energy in the light pulse. It is 7.8 nm for vacuum ultra-violet light (Hu = 100 eV), 
but 1.8 x 10 4 nm for X-rays (Hlj = 5000 eV). We conclude that deep penetration of light 
wrecks havoc with the surface of optical substrates. 



VII. CONCLUDING REMARKS 



We have presented the description of a complete Mesoscopic Hydro-Thermodynamics 
of phonons, which can also be referred to as Higher-Order Nonlinear Generalized Hydrody- 
namics of phonons. It significantly extends the standard hydrodynamics of phonons by 
introducing a complete description in terms of the densities of phonon quasi-particles and 
of the energy, accompanied with their fluxes of all orders [Cf. Eqs.f llOj) and OlSp ]. That 
is, as indicated in the main text, we can talk of the motion of two families, the one as- 
sociated to the motion of the quasi-particle density, together with the evolution equations 
for its fluxes of all orders, and the one for the motion of energy density, together with the 
evolution equations for its fluxes of all orders [Cf. Eqs. ( )29|) and ( 130|) ]. coupled together by 
cross-correlations describing thermo-striction effects. 

As already noticed in the introduction MHT allows to cover all kinds of motion, in that it 
includes those characterized by short wavelengths and ultrafast time evolution. The system 
of coupled equations of motions is extremely cumbersome, in principle of unmanageable pro- 
portions. For handling it is required, depending on each case being considered, to introduce 
a contraction of description implying in retaining only a few number of fluxes, neglecting 
those that become negligible in time intervals smaller than the experimental resolution time. 



In other words 



43], the contraction of description implies in retaining the information con- 



sidered as relevant for the problem in hands, disregarding nonrelevant information. In the 
main text (Section 2) it has been discussed a criterion for deciding on the order of contrac- 
tion. It has been considered in full detail the case of a MHT of order 1, which corresponds 
to a large class of practical situations. The four characteristic Maxwell times are evidenced ( 
they are all important in determining the order of contraction of description). The coupled 
hyperbolic Maxwell-Cattaneo-like equations for the densities of quasiparticles and energy 
are derived, and neglecting thermo-striction effects both sets are decoupled and analyzed, 
obtaining for the heat transport a Guyer-Krumhansl-like equation. 
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Finally in Section 6 thermo-striction effects are taken into account to study the expected 
thermal distortion in silicon mirrors under incidence of high intensity X-ray pulses in Free- 
Electron-Laser facilities, and establishing limiting conditions. 
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Appendix A: The Non-equilibrium Statistical Operator 



According to NESEF ( 12j, 14[, 19] with a short overview given in 15]), the nonequilib- 
rium statistical operator in terms of the basic nonequilibrium variables in sets ([5]) and (J7J) 
is given by 

(t) = Qe (t) X Q B , (Al) 



where 



g e (t) = exp {lng(t, 0) - f dt' e e ^ 

I J — oo 

4ln 5(t ', t '- t )} 



(A2) 



with q (t, 0) being the auxiliary statistical operator (also called "instantaneous quasi- 
equilibrium operator") and 

p(t', t'-t) = expi-m - £ q [> q (0 v^t' -t) + Eq^qQ (0 ? qQ (f " *)] } ( A3 ) 

where t' is the dependence on time of the non-equilibrium thermodynamic variables F and 
the dynamic microvariables, in the Heisenberg representation, depend on (t' — t). Moreover, 
qb is the canonical distribution of the bath of acoustic phonons in equilibrium at a temper- 
ature To, and (j)(t) ensuring the normalization condition plays the role of the logarithm of a 
nonequilibrium partition function. 

We recall that the second term in the exponent in Eq. ( ]A2j) accounts for historicity 
and irreversibility in the non-equilibrium state of the system. The quantity e is a positive 
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infinitesimal that goes to zero after the trace operation in the calculation of averages has 
been performed. We also recall that 



Qe(t) = Q(t,0) + Q' £ (t) (A4) 

i.e., it has an additive composition property, with a contribution of the instantaneous quasi- 
equilibrium statistical operator plus the one of g' E (t) which contains the historicity and 
produces irreversible evolution. 

Next, consider a change of description consisting into going from the one in terms 
of the set (jSJ) to one in terms of the (hydrodynamic-in-character) basic microdynamical 
variables 

[n(r,t), I„(r,t),...,^(r,t),...,Mr,t), I fc (r, t), f \v, t), ...} , (A5) 

whose average values, taken over the non-equilibrium ensemble, of these operators, are the 
macrovariables in sets ( ITUl) and ( 1T51) . which we have dubbed as the n-family and the /i-family 
respectively. In order to calculate the averages we use Eq.(A.3) once we write 

F q (r, t) = ifn (r, t) + F n (r, t) • V q u; q + £^ 2 Fjf (r, t) V^ q 

+ ip h (r, t) tku^ + F h (r, t) -^ q V q ^ q + Ee^h ( r . *) © ^qV^w q , (A6) 

where V q a; q is the group velocity of phonons with crystal momentum q and V q ^ q is given 
in Eq. ffT4l) . The symbol stands for fully contracted product of tensors. 

"Contraction of description" in a given order, say k, is done by taking as null the 
quantities (r, t) and (r, t) for all £ > k. In the main text we have introduced a 
study of MHT of order 1, i. e., keeping only ip n (r, t), (ft (r, t), F n (r, t) and F^ (r, t), and the 



closure of the evolution equations was done using the Heims-Jaynes 44] method as described 
in Appendix C. 

Appendix B: Generalized Peierls-Boltzmann Equation 



As indicated in Eq. (122ft . the evolution equation for the single particle distribution 
^ qQ (t) = (i/i) _1 Tr{[P q Q, H]q £ (t)}. (Bl) 
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Applying the NESEF-based kinetic theory we find 

du^(t)/dt = J™(t) + jW(t) + J®{t), 



where, 



J$(t) = (<«)- 1 rr{[P qQl ff O 5]ff(t,0)} l 
J$i(t) = (iHy'Triiu^, (HsB + Hsp)]Q(t,0)}, 



(B2) 

(B3) 
(B4) 



J^(t) = (my 



(B5) 



Tr{[[H SB + Hsp) , p qQ , [Hsb (0 + (f)JJ]fi (*» 0)} 

Here, 5 stands for functional derivative. We stress that this expression corresponds to an 
approximation where the interactions Hsb and H$p are retained only up to second order 
(memory and vertex renormalization effects are neglected). After performing the calculations 
it follows that 

dVq(}(t)/dt = i (Wq+Q/2 ~ W q _Q/ 2 ) ^qCj(t) + * ( n q+Q/2 ~ n q -Q/2) fqCj(*) 

- (i/2) (r q+Q/2 + r q _ Q/2 ) !/ qQ (t) + j^\t), (B6) 

for Q ^ 0, and 

dv^tyat = -r q (z/ q (t) - u°) + j^\t), (B7) 

for Q = 0, where is® is the distribution in equilibrium at temperature T . Moreover T q , the 
reciprocal relaxation time of the phonons in mode q, is given in Eq. (l25|) in the main text, 
and n q , the self-energy correction of the energy o; q , is given in Eq. fTSH]) . 

For the sake of simplicity, we introduce an approximation of long wavelength, i.e., 
we consider Q « Qb, Qb being the Brillouin radius, in the form 

(B8) 

(n q+Q/2 - n q _ Q/2 ) w Q-v q n q , (B9) 
(r q+Q/2 + r q _ Q/2 ) « 2r q . (bio) 

Then, when going over to direct space, through the Fourier transform of variable Q 
into r, there follows Eq. (l23|) which has a form resembling the standard Peierls-Boltzmann 
equation. 



(^q+Q/2 — W q _Q/ 2 ) ~ Q-V q 0; q , 
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Appendix C: Application of Heims-Jaynes Formalism 

Writing for the auxiliary statistical operator (cf. Eq. (lA3j) 

g(t,0) = exp|l+5}/Tr|exp|l+s}}, (CI) 

where 

A = £ q F q (*) K (C2) 
5 = E q E Q ^q(Q,*)A,Q, (C3) 

that is, a separation in terms of the homogeneous contribution, A, and the departure from 
it, B, and introducing 

g(t,0) = exp {1} /Tr {exp jl}} , (C4) 

the ho mog eneous contribution, according to Heims-Jaynes perturbative expansion for av- 
erages 44] keeping only the first order contribution in the inhomogeneous part of B, and 
using the contraction in MHT of order one, i. e. [see Eq.( j37l) ] 

F q (Q, t) = tp n (Q, t) + ip h (Q, t) ftw q + F n (Q, t) ■ V q w q + F h (Q, t) -/zo; q V q u; q (C5) 

it follows that 

'W*) = [ffl (Q> *) + fh (Q, t) HUq + F n (Q, t) "VqWq + F h (Q, t) -frWqVqWq] 

x i/ q (t) [1 + i/ q (t)] , (C6) 

where 

z/ q (t) = [exp{F q (t)}-l]- 1 , (C7) 

with 

(t) = <^n (*) + ^ (t) + F n (i) -V q O; q + F h (t) -^VqWq, (C8) 

and we notice that 

i/ q (t) [1 + i/ q (t)] = Su^tySF^ (t) . (C9) 

Moreover, we write 

i/ q (t) ~ F q (t) - F q (t) [1 + I7 q (t)] [F n (t) -V q w q + F ft (t) -fku^^] , (CIO) 

where 

F q (t) = [exp {ip n (t) + <p h (t) fko^} - I}' 1 , (Cll) 
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that is, a first order Taylor expansion in F n and (linear approximation) . 

Next, resorting to the use of the nonequilibrium equations of state that relate the four 
nonequilibrium thermodynamic variables in Eq. (IC6p to the four basic variables, i. e., in 
reciprocal space 

n (Q> *) = E^ Q (t) = A u (t) (Q, t) + A12 (t) <p h (Q, t) , (C12) 

I n (Q, f) = 4a (*) • F„ (Q, t) + 41 (*) • Fh (Q, t) , (C13) 

h (Q, t) = An (*) <p n (Q, t) + A 22 (t) <p h (Q, t) , (C14) 

I h (Q, t) = 41 (*) • F„ (Q, t) + 41 (*) " (Q, t) , (C15) 

where A n , A 12 , A®, lf 4 , A 12 , A 22 , A® and are those of Eqs.(|D5]), flDgJ), (JD8]), (JD9]), 



f lD5|) . f lD7j) . f lD9|) and flDTUj) in Appendix D, except for the replacement of z/ q (t) of Eq. flCTOj) 
by F q (t) of Eq.flCTll. 

In Eqs. flC12p and (I014[) the contributions in F n and F h present in Eq. flClOp are null, 
whereas in Eqs. (IC13"|) and (1C15P are null the contributions in ip n and <fh- Eqs. (lC12p to 
( 1C15P constitute a set of linear algebraic equations that can be inverted to obtain the four 
nonequilibrium thermodynamic variables (p n , (ph, F ra and F/j, in terms of the basic hydrody- 
namic quantities, n, h, I n and 1^. 

The second-order fluxes are given by 

J? (Q, t) =4 2 3 (*) (Q, t) + 41 (t) ^ (Q, t) , (C16) 



(Q, t) =A& (t) Vn (Q, t) + A£ (t) <p h (Q, t) , (C17) 

j2] [2] [2] 

where A 33 , A u and A u are those of Eqs. (lD8p . f]D9|) and (IDlOp in Appendix D, except for 
the replacement of i/ q (t) of Eq. flCTOj) by F q (t) of Eq. flCTll . 

Using the nonequilibrium equations of state, after going over direct space we arrive at 
the expressions for the divergence of both second-order fluxes given in Eqs. (H9|) and (|50|) . 
and then to the closed system of Eqs. fl53l to f}56l) . 

On the other hand, introducing the concept of nonequilibrium temperature, better called 
quasitemperature T* (r, t) in the form 

k B T* (r,t) = cp^ 1 {r,t) , (C18) 
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we can obtain an evolution equation for it starting with the evolution equation for the energy 
in the form of the hyperbolic Maxwell-Cattaneo equation, Eq.f J75|) . from which together with 
the nonequilibrium thermodynamic equation of state, Eq. (lC14p . we have that 



E q (^q)XW [l+* q (t)] 



gVfe (r, t) 
dt 2 



+ (Ok 1 + « 



n dip h (r,t) 



dt 



+ V 



■ Vy? h (r,t) 



4 2 i (*) 



A u (t) 

and, after introducing the heat capacity 



A n (t) 

V«(r,t) + rf xt -(r,t) 



(C19) 



Cv (t) = (|^) 2 * q (t) [1 + VM , (C20) 

where T is the temperature in equilibrium in this linear treatment, we arrive at Eq.( l95|) . 



Appendix D: NESEF Kinetic Theory and Expressions of the Coefficients in Eqs. (l42|) 
and (|15D to ([56]) 



The several coefficients for Vqct(t) in Eq.((42]) follows: 

61 (q, t) = A^ 1 (t) [A 22 (t) - hw^A^ (t)] z/ q (t) [1 + z/ q (t)] , 
6 2 (q, t) = A^ 1 (t) [A u (t) fiu; q - A 12 (t)] i/ q (t) [1 + z/ q (t)] , 

b 3 (q, = Agt (t) [41 (*) - ^q41 (*)] ' (*) [1 + ^q (<)] V q W q , 
b 4 (q, t) = A3-; (t) [41 (*) ^q - 41 (*)] ■ ^q (*) [1 + ^q (*)] V q C^ q , 

where the quantities Aij and are, 

^11 W = Eqfq(*) [l+^q(*)], 
A 12 (t) = A 21 (t) = £ q Z/ q (t) [1 + I/ q (t)] fiW q , 



^22 (t)=Eq"q(*) t 1 + ^q (*)] 



q/ ' 



41 (*) = E q ^q (*) [1 + ^ (*)] [V q w q V q a; q ] , 



(Dl) 

(D2) 
(D3) 
(D4) 

(D5) 

(D6) 
(D7) 
(D8) 
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41 (0 = 41 (*) = E q ^ (t) [1 + ^ (*)] [V q o; q V q a; q ] tkv v (D9) 

41 (*) = E q ^q (0 [1 + ^ (*)] [V q ^ q V q o; q ] (^ q ) 2 , (D10) 

A12 (t) = A u (t) A 22 (t) - A 12 (t) A 12 (t) (Dll) 

a 34 (o = 41 (*) © 41 (*) - 41 tt) © 41 (*) • (Di2) 

In these expressions, [V q w q V q w q ] denotes the second order tensor with components 
dco/dqidu/dqj, while © = J2ij F ij G jh and 

z/ q (t) = [exp{F q (t)}-l]- 1 (D13) 

is the population in mode q (see Appendix C), and finally 

*S (*) = E q [bs(q,0 V q n q ] , (D14) 

41W = E q [b 4 (q,t)V q n q ], (D15) 

41 (t) = E q ^q [b 3 (q, *) V q II q ] , (D16) 

41 (0 = E q ^q [b 4 (q, *) V q n q ] , (D17) 

6n(*) = -E q 6i(q,*)r q , (Dis) 

6i2(*) = -E q 62(q,*)r q , (Dig) 

M0 = -E q &i(q,0r q ^ q , (D20) 

6 22 (0 = -E q &2(q,0r q ^ q , (D21) 

41 (*) = "E q [V q a; q b 3 (q, f)] T q , (D22) 

41 (*) = -E q [V qWq b 4 (q, t)} T q , (D23) 

41 (*) = "E q [V q ^ q b 3 (q, *)] T q ^ q , (D24) 

41 (*) = "E q [V q o; q b 4 (q, *)] T q ^ q , (D25) 

and we recall that we are writing [AB] for the tensorial product of vectors A and B, 
rendering a tensor of order two. 
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